NASA Contractor Report 189665 
ICASE Report No. 92-24 





A COMPUTATIONAL ALGORITHM FOR CRACK 
DETERMINATION. THE MULTIPLE CRACK CASE 


Kurt Bryan 
Michael Vogelius 


Contract Nos. NAS1-18605, NAS1-19480 
June 1992 

Institute for Computer Applications in Science and Engineering 
NASA Langley Research Center 
Hampton, Virginia 23665-5225 

Operated by the Universities Space Research Association 


IVIASA 

National Aeronautics and 
Space Administration 

Langley Research Center 

Hampton, Virginia 23665-5225 


vC 

f-4 

UO 


sO 

O 

in 

jTi 

m 

m 

o 

1 




O 


a* 

c 

fH 


’"'I 

c 


& 

nO 

\ 

m 

o 


< e 

C H* 
*-* <t 


<X v — ■ * 


IT 

C 


D a lL 
Q U i 
X *- oj 

C' U‘ 00 

u a < a 
cj 

< — * 

lj y 
< o 

^ o: <r 

LCt oj or ^ 

<0 Li . 

O ur - O 

C’ f u < 


r~* a ►- 

t >: * * ^ 

0 rr t — 
j t- 

1 D *-• 
-t y u 

i 1 O 

< o .u 
— IT Cl 
^ < I— .L 



A COMPUTATIONAL ALGORITHM FOR CRACK DETERMINATION. 


THE MULTIPLE CRACK CASE 1 


Kurt Bryan 

Institute for Computer Applications in Science and Engineering 
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Michael Vogel i us 
Department of Mathematics 
Rutgers University 
New Brunswick, N.J. 08903 


Abstract 

This paper develops an algorithm for recovering a collection of linear cracks in a homogeneous 

electrical conductor from boundary measurements of voltages induced by specified current 

fluxes. The technique is a variation of Newton’s method and is based on taking weighted 

averages of the boundary data. The method also adaptively changes the applied current flux 

at each iteration to maintain maximum sensitivity to the estimated locations of the cracks. 

l The first author is partially supported by National Aeronautics and Space Administration contracts 
NAS1-18605 and NAS1-19480; the second author is partially supported by National Science Foundation 
Grant DMS-89-02532 and AFOSR Contract 89-NM-605. 
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1 Introduction 


In this paper we develop a very efficient computational algorithm to reconstruct a collection 
of linear cracks inside a homogeneous conductor from electrostatic boundary measurements. 
The algorithm in this paper can be seen as a natural extension of the algorithm developed in 
[15] for the reconstruction of a single crack. This extension poses several theoretical as well 
as practical challenges. We have also significantly improved the efficiency and versatility of 
the earlier algorithm by basing all computations for the underlying conductance problem on 
a one-dimensional boundary integral formulation instead of a two-dimensional finite element 
formulation. It should be mentioned here that boundary integral formulations have been used 
in other implementations of (single) crack reconstruction algorithms [13], [14], and also that 
progress on the development of an algorithm for the reconstruction of a single (penny-shaped) 
crack inside a three dimensional object is reported in [14]. Algorithms like the one discussed 
in this paper are significantly different from more general purpose imaging algorithms (c/. 
[3], [6], [9], [16]) that seek to reconstruct an unknown distributed conductivity profile from 
similar boundary measurements. Algorithms such as ours are based on the assumption that 
certain apriori information about the profile is available and they incorporate this knowledge 
into the reconstruction in such a way as to achieve better continuous dependence and better 
discrete approximation properties. One important feature of the present algorithm is that 
it is based on an adaptive change of the prescribed boundary current patterns to ensure 
“maximal” sensitivity. The idea to use some kind of “optimal” current pattern in connection 
with impedance imaging has been developed by Gisser, Isaacson and Newell (c/. [9]); the 
specific strategy we use is somewhat different from theirs and ties in directly with the iterative 
procedure (it does not rely on any eigenfunctions). Our reconstruction is based on the usage 
of relatively few averages of the boundary voltage measurements (as opposed to all the 
boundary voltage data). In addition to improving the efficiency of the algorithm, this should 
also decrease the probability of getting caught in a local minimum when compared to more 
standard output least-squares algorithms. 

An outline of this paper is as follows. In Section 2 we present the “customary” mathe- 
matical model of electrostatic conductance for a smooth, isotropic background medium that 
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contains a collection of cracks. In particular we demonstrate the duality between the notions 
of perfectly insulating cracks and perfectly conducting cracks. We also briefly discuss known 
uniqueness and continuous dependence results. Section 3 contains a detailed description of 
the 4n functionals that we use for the reconstruction of a collection of n or fewer linear cracks. 
This section also provides a discussion of the adaptive strategy that we use for the selection 
of the “maximally sensitive” electrode locations. The boundary integral formulation of the 
electrostatic conductance problem and its discretization by Nystrom’s method is the topic of 
Section 4. The central part of the reconstruction algorithm is a version of Newton’s method. 
This particular version together with the required gradient computation is discussed in Sec- 

I 

tion 5. Section 6 contains a selection of representative computational experiments with our 
algorithm. Finally in Section 7 we provide a brief summary of our results and a description 
of possible future developments. 


2 The Mathematical Model 

A single crack is commonly modeled as a perfectly insulating curve a. With a background 
conductivity 0 < 70 < 'y(x) < 71 and a finite collection of cracks £ = U£ = 1 < 7 *;, the steady 
state conductance equations thus read 

n\£, (2.1) 

v = <f> on <9fb ( 2 . 2 ) 

The field v is normal to E. The function v represents the voltage potential induced in fl. 

We assume that 0 is simply connected, be., it has no holes, and so the entire boundary 30, 

is accessible from the “outside”. Let u denote the “ 7 -harmonic” conjugate to v. It is related 

to v by the formula 

(Vu) v * * * * x = 7 VU, (2.3) 


V • ( 7 V 1 ;) = 0 in 
dv 

i- = 0 ° n 


with appropriate boundary conditions on dft, e.g., 
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where _L indicates counter-clockwise rotation by 7 t/ 2. For a particular set of constants 
Cfc, k — 1, ... , n, the function u solves the problem 

V-( 7 _1 Vu) = 0 in fi\E, (2-4) 

u = c k on a k , k = 1, ... , n 

with 

7 -1 ^ = t/> = ^ on dfl. (2.5) 

(71/ C7S 

Here s denotes the counter-clockwise tangent direction on $0 and v denotes the outward 
normal on dfl. For these particular constants, finding a solution to (2.1), (2.2) is thus 
equivalent to finding a solution to (2.4), (2.5). The constants c k may (up to a common 
additive constant) be characterized in several equivalent ways: 

(a) Let A be the n x n-matrix with elements Aij = f u ~y~ l X7U^S7U^dx and let b be the 
n-vec.tor with elements 6; = f^QipU^ds, where U^\ j — 1 ,. . . ,n, denote the solutions 
to 


V- ( 7 -'W (j) ) = 0 in ft\£, 

U U) = 1 on er,, 

f/bl = 0 on a k , k ^ j 

with 

7 1 — - — = 0 on dll. 
dv 

Then the vector c — (cj, .... , c n ) ( is the solution to Ac = 6. 

(b) The “ 7 -harmonic” conjugate, u, satisfies / (Tj .[ 7 -1 |^] ds = 0, 1 < k < n. Here 

= 7 -i Jjl _ 7 -1 ^- denotes the jump in the normal flux across the curve a k . 2 
Furthermore, the set of constants { Cfc }*=i * s unique se f constants for which the 
solution to (2.4), (2.5) has this property. 

2 The expression denotes the limit of the derivative (in the direction v) as one approaches a k from 
the side to which e points. denotes the limit as one approaches a k from the opposite side. 
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(c) Let T be a fixed point on < 9 f i, in a neighborhood of which <f> is smooth. Let r k be 
a smooth curve in 0 \ E connecting T to an interior point of the crack <7*., and let 
s denote the unit tangent direction along 77, pointing from T towards 07. Then the 
constants c k are given by the formulae 

ci ‘ = -i k ' i £ ds+u(n 

where u denotes the normal field v = — s 1 . 

The characterization (c) is a direct consequence of the relation ( 2 . 3 ). The characteri- 
zations (a) and (b) are practically much more useful; they are both a consequence of the 
following well known result from convex duality. 

Proposition 2.1 If <f> is an element of then the field £ = -/Vv is the (unique) 

minimizer of the functional 

0 / 7 _1 M 2 dx - [ <j)T]-uds (2.6) 

1 Jn Jan 

in the space H = L 2 {Vt) D {77 : V • t] — 0 in \ £ , 7 • v - 0 on 07, k=l v . . ,71}. 

It is not difficult to see that any element of the space H satisfies V * rj = 0 in all of 
f 2 and therefore has the form 7 / = (Vic ) 1 for some w £ H l ( fi), with w being constant on 
each 07. Conversely, it is also true that any vector field of the form 7/ = (Vic) 1 , w 6 H l (fl) 
n{ic = constant on each 07, k = 1 , ... , 7i}, is in the space H. After insertion into ( 2 . 6 ) we 
obtain 7VU = (Vu) 1 , where u is the minimizer of the functional 

— f 7 _1 |Vic| 2 dx — ! ds 

l Jn Jdn os 

in the space // 1 (fi) n{tu = constant on each 07, k = 1, ... , n}. This provides a variational 
characterization of the ^-harmonic” conjugate to v. Let F(d), d 6 lR n , denote the expression 

m = U r'\^ + t , d ^ ,k) )\ 2 d x - / 

z JQ k = 1 Jau Jt=l 

From the above discussion it is clear that the set of constants corresponding to u are char- 
acterized by the fact that d = 0 is a minimum of F. Equivalently (because of the form of 
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F), the set of constants corresponding to u are characterized by the fact that d — 0 is a 
stationary point for F. Stationarity of d — 0 is equivalent to the conditions 

/ 7 - l S7uVU {k) dx = [ ipU {k) ds , k = 1 , ... ,n. (2.7) 

Ja JdQ 

From integration by parts (on the domain ft \ S) we have 

/ 7 _ 1 VuW (fc) dx = / [ 7 -1 7T~] ds+ / 
in ' icr* 01 / ian 

Insertion of this into (2.7) now gives that the set of constants corresponding to u are char- 
acterized by 

/ ds = 0 k = 1 , ... ,n, 

as asserted in (b). 

On the other hand, the function u has the form 

u = u 0 + £c i U", ( 2 - 8 ) 

i = 1 

where no is the solution to 

V • ( 7 -1 Vu 0 ) = 0 in 
u 0 = 0 on 

with 

7 -i ^2 = 0 on dkl. ( 2 . 10 ) 

ov 

Because of (2.9) and (2.10) we have f n 7 _1 Vn 0 Vf/ (fc) dx = 0; by insertion of (2.8) into (2.7) 
and use of this formula we now obtain 

Yd I 7 ~ 1 Vt/ ( ‘ ) Vf/ (fc) dx = f ipU {k) ds k = 1 , . . . , n, 

“ in Ja n 

which is exactly the characterization (a). The above argument rests on the fact that <f> is 
in H x ^ 2 (dkl) (ip is in H~ x ^ 2 (dkl)). However, by continuity the characterizations carry over 
to cases in which <f> (and ip) are not necessarily so regular that u is a variational solution. 
In particular these characterizations remain valid if ip consists of delta functions, a type of 
boundary current we shall repeatedly use in this paper. 


ft\E, 

k 1 , . . . , tj, 
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Since u and v are related by the equation (Vu) x = 7VU, it is clear that knowledge of the 
pair ls equivalent to knowledge of the pair It is much more convenient 

to work with the function u as opposed to the function v ) and we shall entirely do so for 
the development of our algorithm. In particular by working with u we avoid the difficulties 
that are associated with non integrable kernels in the integral equation formulation ( cf ’ [12], 

[ 14 ])- 

Let Po, ... , Pm be M + 1 points on dfl; we assume that these points are labeled in order 
of counter-clockwise appearance, starting from P 0 . For the crack reconstruction we utilize 
solutions corresponding to the two-electrode currents ifj = 6 Po - 6 P} , j = l,...,Af, 

V • (7 -1 Vuj) = 0 in fl\E, ( 2 - 11 ) 

Uj = cj^ on crk, k = 1, ... ,n, 

with 

7 1 ^ on 9Q, (2-12) 

the constants being selected so that /<,*[ 7 -1 ^] ds = 0 for A; = 1, ... ,n. The inverse 
problem may now be stated explicitly as follows: 

We seek to reconstruct the collection of cracks S = from knowledge of the boundai'y 

voltage data corresponding to the prescribed two-electrode currents 'y - 1 -7^- = 6p 0 ~ dp., 

J = 1 , ... ,M. 

It is known that boundary voltage measurements corresponding to M = n T 1 fixed two- 
electrode currents suffice to uniquely identify a collection of n (or fewer) cracks [4]. This 
result is an extension of a result in [8] which asserts that boundary voltage measurements 
corresponding to two fixed two-electrode currents suffice to uniquely identify a single crack. 
Recently an interesting continuous dependence estimate has been obtained for the case when 
the background conductivity is constant and there is at most one crack [1], Briefly described, 
this estimate states that if the boundary voltage data (on some open subset of dfl) deviate 
by e then the crack locations differ by at most [log(| log e|)] -1 / 4 . In the present paper shall 
always try to fit the data entirely by means of linear cracks. For such cracks one can 
hope to have better continuous dependence estimates, as indicated by the results in [8]. As 
was the case in [ 15 ], we base the reconstruction of the cracks on the values of a number 
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of functionals (as opposed to all the boundary voltage measurements). In [15] we used 4 
functionals corresponding to the reconstruction of a single linear crack; the natural extension 
is to use 4n functionals for the reconstruction of n cracks. In the following section we give 
a careful description of these functionals. 


3 The Functionals 

We now specialize to the case of a constant background conductivity, 7=1. Let F denote 
the vector- valued function 

F(E,V>,w) = (F(E,V , ,t« (1) ),F(E,V>,^ (2) ), F(E,V», w (3) ), C(E, 0, w (4) )) ( , 
where F(E,i/>,u;) is given by 

C(E, w) = u(E, 4’)-^ ds, (3.1) 

and where 1 < i < 4, are particular solutions of 

Aw = 0 in 1R 2 \ £. 


The function u = u( S, 0) is the solution to 


A u — 0 in fi\E, 


u 

du 

du 


Ck on cr*, k — 1, . . . , n, 
ip on 


(3.2) 


with the constants c* uniquely specified by 


/ u ds = 0 and / 
^30 Jerk 


du 


du 


ds = 0, A; = 1, . . . , n. 


The exact selection of boundary currents and test functions w = (iy^\ u/ 2 ), w^ 3 \ u;^)* is 
very important and will be discussed shortly. We select one tp and one w corresponding to 
each crack er*; whenever we want to emphasize this correspondence we use the notation 


xp 




and 


w 


°k 


— ( w , 


(!) w m 

°k » W <*k 


W 


( 3 ) 
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We have for convenience chosen the normalization f da u ds — 0 for the voltage potential. We 
shall always select w so that 



dw ^ 
dv 


ds = 0, 


and 



ds — 0 k — 1, . . . , n. 


(3.3) 


Because of the first identity in (3.3) the function F is unchanged by the addition of a 
constant to u , and we could therefore just as easily work with any other normalization. 
The components of F are just weighted averages of the boundary voltage data. We use a 
weighting function of the form —j because of the relation 

f dw f • 

/ u(E,0)^r— ds = / Vu(E, xp)Vw dx (3.4) 

JdU dv J n\E 


that exists between the expression in equation (3.1) and the energy bilinear form (by means 
of Green’s formula). As will be seen later, these averages are equivalent, in the absence of 
any crack, to the set of the first 4 nontrivial Fourier modes of the induced boundary voltage. 

The data for our reconstruction consist of measured boundary potentials corresponding 
to certain prescribed two-electrode boundary currents. We denote the voltage data cor- 
responding to the boundary current xp by g(xp) and define a corresponding vector- valued 
function 

f (V>, w) = (/(0,u> (1) ), /(</>, w (2) ),/(^, t« (3) ),/(V’,^ (4) )) < , 


where f(xp,w) is given by 

f{ip,w) = f g ds, 

JdU OV 

and where are the same functions as before. Our algorithm seeks a solution E = {cr^} ? fc l =1 
to the 4n equations 

F(E,Vv fc , w,J = f (0<r*,w„J, 1 <k<n. 


Consequently, we do not use information about the full boundary voltages for the recon- 
struction; we only use information about the values of these particular functionals. 

We implicitly assume that our data is consistent so that f (xp ( Tfc , w^ fc ) corresponds to some 
collection of cracks E* (E* may consist of fewer than n cracks). In reality we therefore solve 


Gfc(E) = 0, 1 < k < n, 


(3.5) 
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with 


Gfc(E) = F(E 1 ^,*J-F(E* l K,w, l ) (3-6) 

= F(E, w^) -f (K.WaJ- 

Clearly E* is a solution to (3.5). If the Frechet derivative D E {Gfc}|s = £* (a 4nx4n matrix) 
is nonsingular, then E* is indeed the unique solution to (3.5) near E*, and furthermore, one 
may expect that some variation of Newton’s method will be an efficient solution technique. 
Differentiation with application of the “chain rule” yields the expression 


Dz{G k }l = :1 |e = e* = {£> e F(E,i/v*, w„ ; )| s= £.})E =1 , 

for the derivative with respect to E (at E*), the right hand side of which is a 4n x An matrix 
with rows 


dw 


(0 


D e F(E,^*,</)|s=e* = / 8n £>E«(S,^;)|E=E-5^ 

^ K 

1 < k < n, 1 < f < 4. 


ds, 


(3.7) 


In [15] we explicitly calculated the expression (3.7) in terms of u and (eliminating D E u), 
we used this alternate expression for the selection of the “maximally sensitive two-electrode 
currents as well as for the Newton’s update itself. In our implementation here we shall rely 
on essentially the same technique as in [15] for the selection of two-electrode currents, but 
we shall directly compute the derivative of t/ ( E , ) with respect to E at the discrete level 

(cf. Section 5) for use with the Newton’s update. 

When talking about the derivative with respect to E, we mean the derivative with respect 
to the 4n parameters that are used to describe E. Just as in [15] we parametrize a single 
crack by (6,, b 2 ), 0 , and A, where (b u b 2 ) are the coordinates of one endpoint, 6 is the counter- 
clockwise angle between the crack and the halfline y — b 2 , x > £>i , and A is the length of the 
crack (here coordinates of points are relative to a fixed reference coordinate system). For 
convenience we order these coordinates q = (0, b 2 , &i, A); the parameter set corresponding to 
E is now given by the vector Q = (qi,q 2 , ... , q») ( - The derivative D^{G k ] consists of the 
(n 2 ) 4x4 blocks 

D qi G k . 
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At the solution, E = £*, these blocks equal 


Z) qi F(Q,^.,w CT .)| Q=Q .. (3.8) 

If the derivatives (3.8) are formed at some Q° not equal to Q* then they no longer represent 
the full derivative of {G*}; the latter also includes terms where w (7k and ift (Tk are differentiated 
through (c/. (3.6)). For the Newton’s update we use a more complete “derivative” which 
includes the terms that arise when the are differentiated through, but we do not include 
differentiation through t/> ak . The goal behind the choice of w crk and ip (7k is to make the 
derivative of {G fc } as far from singular as possible at E =f E°, the current stage of the 
iteration. 

In order to describe the choice of we select a coordinate system such that er *. lies on 
the positive x\ axis, with one endpoint at the origin. In this coordinate system we choose 

= Im ( 2 L w ll ] = Im[z 2 ], (3.9) 

Re[(z - A k)Jz(z - A*)], Re(z ) > 

, (3-10) 

-Re[(z - A k)\jz{z - A fc )], Re(z) < ^ 

Re[v /^( 2 ~ A*,)], Re[z) > 4^ 

, (3-11) 

-R e[yjz{z - A*)], Re(z) < ^ 

where z — x\ -\-ix 2 and A*, denotes the length of cr*. The functions and are extended 
to Re(z) = Afc/2 by continuity. Except for a change in this is exactly the same choice 
of test functions as in [15]. Since they are harmonic in fl, the two functions and wjfj 
clearly satisfy (3.3). It requires some extra calculations to check that wf^J and w^ 4 J also 
satisfy (3.3); for reasons of brevity we omit these calculations here. The fact that uA 3 ) and 
w ^k > n ^ ee< 3 satisfy (3.3) makes the Remark 1 in Section 3 of [15] superfluous. Notice that 
w i^> 1 5; i < 4, vanish on the crack cr*. 

Remark 

Given the specific form of the weight functions uA*] it is now fairly easy to explain why, in 
the definition of G*,, we pick a distinct boundary current xj) ak corresponding to each k. If as 
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an extreme we had picked the same boundary current rf corresponding to each crack, then 
the first two equations of (3.5), (3.6) would be solved for all k if 

L&MV* = Li {x)sWis ' 
is = Ll {y)9Wds ' 

[ ^-(x 2 - ?/ 2 )u(S,V>) ds = l ^-(x 2 - y 2 )g{d>) ds 

Jan ov Jdn Ov 

r\ f C) 

I — (zy)u(E,V>) ds = I —(xy)g{ip)ds, 

JdU OV Jdtl OV 

where (x,y) denote coordinates relative to some fixed coordinate system. The system (3.5), 
(3.6) would therefore represent no more than 2 n + 4 equations for the An unknowns of E. 
For n > 2 this immediately leads to “underdetermination” and a singular Jacobian. ■ 

In [15] we analyzed the structure of the derivative 

D q F(q,Vv>, w ff o)| q=q o, (3.12) 

for the case of a single crack (and with a slightly different choice of u;^). With the ordering 
of the parameters (0, b 2 , &i, A) we found that this 4x4 matrix was lower triangular. We also 
found that if the test functions w^'\ 1 < * < 4, had been selected harmonic (and 0 on the 
crack) then the last two columns in this matrix would have been zero. Test functions with 
singularities like and ud 4 ^ are thus essential to insure that this matrix is non-singular. 
Since the only change we make in the test functions concern u/ 4 I we do not destroy the lower 
triangular structure of this derivative for the one crack case. In the multiple crack case the 
counterparts of the matrix just discussed are the diagonal entries 


^q*F(Q,^ <r o,W )T o)|Q = Qo. 


It is worthwhile noticing that these matrices do not inherit the lower triangular structure. 
For a any fixed crack Ok contributions corresponding to the other cracks will appear above 
the diagonal. These contributions are of the form — J (T( d s i where the functions j 

are related to w{^ by the formulae 

& k J 


/W*) 

^ 2 , = (x,,x 2 )‘ - Vwi‘1, 
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c/. [15], page 921. However, it is our practical experience that the above selection of 
together with the appropriate selection of the two-electrode current ip Ck (to be discussed 
below) is a very effective way of achieving a Frechet derivative which is far from singular. 

Remark 

It is interesting to consider the limit A*, — > 0 as one endpoint and the direction of the crack 
stay fixed. Let a k be a crack with endpoint at (b,,b 2 ), zero angle, and length \ k . The 
corresponding limits of the functions are given in absolute coordinates by 

w o } = y~b 2, ?4 2) = 2(x - b,)(y - b 2 ), (3.14) 

4 3) = (x - b,) 2 - (y - 6 2 ) 2 , ?4 4> = X - bi. 

Actually, = Wq ^ and w^ 2 J = w t identically, since these two functions do not depend 
on the crack length. Moreover, the approximations Wo 3) and w£> ~ Wq 4 ^ for the third 

and fourth functionals are quite accurate sufficiently far away (e.g. two crack lengths) from 
the crack. If fi is the unit ball and g(9), 0 < 9 < 2ir, is the limiting boundary voltage, and 
Ff = limA*_o Vb «>£*)> then (3.14) gives 

F\ = ft g{Q) sin 0 d9 

F° = 2 ft g(9) sin 29 d9 - 2b, ft g(6 ) sin 6 d6 - 2b 2 ft g (9) cos 9 d9 
F° = 2 ft g{9) cos 29 d9 - 2b, ft g(9) cos 9 d9 + 2 b 2 ft g(9) sin 9 d9 

F! = ft g(9) cos 9 d9. 

In the limit — > 0, F thus represents the first 4 Fourier modes of the boundary data. ■ 

The change we have made in can be explained from this remark: with the choice made 
in [15] u/ 3 ) and had the same limit as A — * 0. The boundary integral implementation of 
the algorithm would therefore occasionally attempt to fit the data by just making the cracks 
very small, since it then in effect only would have to satisfy 3 n equations (using 3n unknowns) 
as opposed to satisfying a larger set of 4n equations. We did not encounter this phenomenon 
in any of the single crack experiments performed in [15] since the implementation there 
was based on a 2-d finite element formulation, which effectively put a lower bound on the 
crack length. The new we have introduced here is simply a (correctly scaled) linear 
combination of the and used in [15]. 
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For the selection of the current xj> ak we rely on an appropriate adaptation of the technique 
developed in [15]. In that paper we calculated that the first and second diagonal entries of 
the matrix (3.12) are m and 2m respectively, where m denotes the expression 

dw\ l ^ 
dv 

We then proceeded to select a two-electrode current which made this expression largest 
possible. In the multiple crack case the first diagonal entry of the k’th diagonal block of 
{£> s F(E,Vv fc ,w^)}£ =1 is given by 



'tfnv^ ] ds. 


m k 


= / (“^r 1 




ds — 


r du 
J a i Ov 


XV. 


( 1 ) 

1 


ds 


(3.15) 


where is as defined by (3.13). 

Consider the function £* 6 H l {Sl \ a^) satisfying 


A£ k = 0 in \ <T fc , 


& = 0 on crjt, 

djk _ 

dv 8v 


(3.16) 


on <90. 


Note that does not identically vanish on cr k , so that £ k is different from After 

insertion of ( k into (3.15) we have 


m k 


= L [ u i£ - ds -^ k L, 

= f M m -41’,) + /„ ([^1& - “[|r ]) ds 


du 

Tv 


(3.17) 


= J g n V>(6 - ) ds + S/^jfe jf (6 - «£>.,) 


To obtain the last identity we have used the fact that £*. = 0 on cr* and that 


llfus = -l £* — / ^*- 0 . 

Ja k dv JdQ dv JdQ dv 

Concerning the two terms in the last expression of equation (3.17), it is reasonable to assume 
that the first term JdQ — w ^k, i) d s be ,ar § er ( a ^ l eas f f° r moderate size cracks). 


13 


If we now substitute xp = dp — 8q into the last expression in (3.17) and disregard the second 
term we get 

"H w (fc - w^ :l )(P) - (& - u.^., )(<?). (3.18) 

We now chose tp ai = Sp 0 — 6 Pl so that P 0 maximizes (£] — )(P) and P\ minimizes 

(£i — | j )(Q). For subsequent 2 < k <. n we chose xp ak — 6p 0 — bp k (the same Pq as for Ar = 1 ) 

where Pjt is selected so that the expression (£*, — w ^ kt x)(Po) — (£* — w ^ k ,\){fk) is of maximal 
magnitude, subject to the constraint that the points P 0 , P \ , P 2 , . . . , P n stay well separated. 
We cannot allow any of the P 3 to coincide, for then we would be duplicating boundary 
measurements, potentially leading to a singular Jacobian as described in the remark before. 
The above construction makes the two-electrode currents appear somewhat like the currents 
used for the uniqueness result proven in [4]. For the uniqueness result we needed n + 1 
two-electrode currents with n + 2 distinct electrode locations - the currents prescribed above 
only number n (with n + 1 electrode locations). We expect that this deficiency is more than 
compensated by the fact that the electrode locations change as the iterations proceed. 


4 Integral Equation Formulation and Discretization 

We now proceed to formulate the boundary value problem (3.2) as an integral equation on 
the boundary of the region fl \ E. Let T(x,y) denote the fundamental solution for the two 
dimensional Laplacian given by 

r(*,y) = 2-l°g(k - y|), x,yelR 2 . 

The application of standard potential theory arguments (see [7], Sections 3.B-3.D) shows 
that u, the solution to (3.2), for x € fi \ E, can be represented as 

«(*) = (( u (y)-S-^(x,y)-T(x,y)^(y))dsy-J2 I ?{x,y) 

J ail OP y 

Here^ on dQ denotes the normal outward derivative with respect to the y variable. On 
any of the cracks <r fc , v denotes a unit normal field and = dhi. — denotes the jump 
in the normal flux across the crack. To simplify notation we shall use the notation <p^ for 


du 


(y)ds y . (4.1) 
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the jump 1 1^| across cr fc . The formula (4.1) expresses the value of u at any point in fi\ E in 
terms of the Dirichlet and Neumann data for u on dfl and the jump in across the cracks 
E. This is simply Green’s third identity applied to the region 1) \ E. It is straightforward 
to check that since has at most an r -1 / 2 singularity at the endpoints of any crack and 
since the singularity of T at x = y is only logarithmic, equation (4.1) also holds for x € E. 
Because u is a constant c; on <t/ (and continuous in il) this implies that 

/ ( u (y)J—F( x ,y) - ^( x ,y)Hv)) ds y - it, [ r(x,y)<f> k {y)ds y = c i (4.2) 

JdQ UlSy J(T k 

for x 6 For x £ di an argument similar to that which led to (4.1) leads to the equation 


\u(x)+ j u(y)^-r(x,y)ds y -J2 [ T(x,y)(f> k (y) ds y = [ T(x,y)tf>(y) 
L JdQ C'L'y — \ JdQ 


ds y . (4.3) 


As discussed earlier, if the constants c/ are treated as unknowns they can be determined 
uniquely from 

f <f)\ds — 0, / = (4.4) 

J a i 

with the normalization f du uds = 0. Combining these n + 1 conditions with equations (4.2) 
and (4.3) we arrive at the following system of integral equations 

d 


1 r a n t 

-~u(x)+ u{y) — r(x,y)ds y - Y(x,y)<}> k (y) ds y 

J JdQ UtSy J<T k 

= / r {x,y)%l?(y)ds y , x G Oft 

JdQ 


(4.5) 


[ <‘(y)-£~r(x,y)dsy 

JdQ 0u v 


Y, F{ x ,y)<!>k(y)dsy - c t 

fc=i Jak 

/ r(x,y)ij>(y)ds y , x £ a t 

JdQ 


with the additional constraints 


j 

JdQ 


<f>l ds 


0, l = 1, . . . n, 


uds = 0. 


(4.6) 


(4.7) 


The unknowns here are the value of u on <9f l, <f> k on a k and the constants c;, l = l,...,n. 
Given a set of cracks E and Neumann data ip one can solve equations (4.5), (4.6) and (4.7) 
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to obtain these quantities. The solution to the boundary value problem (3.2) at any point 
in can then be obtained from the representation (4.1). 

One useful fact to note is the following. Let Uo be a harmonic function on with 
Neumann data xp. The same reasoning used to derive the integral equations for u shows that 
Uo satisfies the boundary integral equation 

1 f d f 

-~u 0 (x)+ u 0 {y) — Y{x,y)ds y = Y(x,y)^(y) ds y (4.8) 

2 JdQ UlSy Jd£l 

for x E <9fL A unique solution again requires a normalization such as f d Q u 0 ds — 0. If x E 
then uq(x) can be represented 

“o(i) = / (u 0 (y)-S-Y(x,y) -T(x,y)ip(y))dsy. (4.9) 

JdQ OUy 

Let v denote the difference u — u 0 . Combining equation (4.5) with (4.8) and combining (4.6) 
with equation (4.9), we find that v satisfies 


■\v(x)+ f v(y)-^-r(x,y)ds y ( r (x, y)^ k {y) ds y = 0, x € dQ (4.10) 

^ (s IS y j J (T 

f <9 n ^ r 

/ v (y)-^—r( x ,y)d s y-J2 Y(x,y)<f) k (y)ds y - c t = -u 0 (x), x € cr, 

JdU UlSy ^ & k 


where <f> k denotes the jump in the normal derivative v across cr^. Note that this is the 
same as the jump in 4-u, since uq is smooth in fi. The conditions 


I (j>ids — 0, / = 1, . . . n, 

J < 7 1 

v ds — 0, 

Jan 


(4.11) 


are also still enforced. Equations (4.10) and (4.11) provide a means of directly computing 
the perturbation v = u — u 0 caused by the presence of the cracks. This formulation is more 
advantageous than the original formulation since we will use singular boundary data xp. The 
function v = u — uo, however, is smooth up to d£l and hence avoids any of the problems 
associated with the lack of regularity of ip. Moreover, for the specific two-electrode Neumann 
data (and a domain in the form of a ball, as considered later) we have a closed form solution 
for u 0 . 

Having derived the boundary integral formulation we now briefly discuss how we discretize 
it by means of the so-called Nystrom method. Suppose that the boundary of the region f) 
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is parameterized by z(t) = (z\(t), ^(O)* 0 < t < 1. Let each crack <7^ be parameterized by 
z(<) for k < t < k + 1. In terms of this parameterization, equations (4.10) and (4.11) can be 
written as 


1 r\ n r k + 1 

-~<t>{s) + K(s,t)ct>(t)dt-Y, G(s, dt = 0, sg[0,1) (4.12) 

2 Jo fc=1 Jk 

r\ n fk + 1 

/ K{s,t)(i>{t)dt- Y, / G{s,tW)dt-c = -u 0 (z(s)), s e [/,/+ 1) 

2o fc=1 2/c 


and 


W + l 

/ /(t)|z'(t)l dt = 0, l = 1, . . . , n, 

f <f>(t)\z'(t)\ dt = 0 
Vo 

where A(.s,f) = g|_r( 2 (.s), z(t))\z'{t)\ and G{s,t) = T(z(s), z(t))\z' (t)\. The functions v and 
<j>k have also been replaced by the single function <j) defined on [0, n + 1] by fat) = v(z(t)) 
for t £ [0, 1), fat) = <j>k(z(t) for t € [k,Jc+ 1 ). 

Let t 3 and u>j, j = 1, . . . ,m, denote the nodes and weights of a quadrature rule on [0,1], 
so that 

j 771 

/( 0 dt » 

3 = 1 

for reasonably smooth functions /(£). Nystrom’s method for solving an integral equation of 
the form 

<£($)+ / = /(*) 

Jo 

consists of replacing the integral by a quadrature rule to obtain 


fas) + = Z( 5 )- 

j=i 

Letting s assume the values ij , . . . ,£ m we obtain the m x m linear system 

m 

(j) t ) ] Kijfa — l 1 , . . . , 777. 

7=1 

where fa = <^>(L), I< t j = K(U,tj)uj and /, = /(*,•). The intention is that fa = <£(f,-) ~ faU). 
A complete treatment of Nystrom’s method for second kind Fredholm equations can be found 
in [2], 
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For a first kind integral equation with a smooth kernel, 

Th'<t>= I K(s,t)<f>(t)dt = f(s), 

direct discretization usually leads to a linear system which is very poorly conditioned. This 
stems from the fact that the corresponding inverse operator is unbounded on whatever 
space the equation is posed, typically <7(0, 1) or L\ 0, 1). The singular values of the operator 
T k approach zero rapidly, so that the solution <f> is extremely sensitive with respect to / or to 
noise on the right hand side of the equation. The smoother the kernel, the faster the singular 
values decay and the poorer is the conditioning of the linear, system. The solution of first 
kind integral equations therefore often requires some kind of regularization. As discussed 
in [5], however, if the kernel K(s,t ) is singular enough on the diagonal s = t, reasonable 
results can be obtained from a direct discretization of the equation, without regularization. 
In the case of equations (4.12) the first kind portion of the equations (on the cracks) have 
a logarithmic singularity along the diagonal a = t, and hence regularization should not be 
necessary. We have applied Nystrom’s method directly to the equations. The linear systems 
obtained in this manner have good conditioning and in all test cases in which we have a 

closed form solution, this method has produced solutions in complete agreement with the 
closed form solution. 

To apply Nystrom’s method to the equations (4.12), we replace the integrals over the 
intervals [/, / + 1], / = 0, . . . , n, with the quadrature rule and then let s assume the values 
l + t t , i - 1, . . . ,m, / = 0, . . . , n. This yields the following linear system in the variables 

<^01, • - • , (j>0 m, <j> 11, , <j>nm,Ci , . . . , C n \ 

1 n m 

2^o» + 22 h (ti, tj )utj <f>Qj — '22 *22 G(t t , k + tj)uj(j)kj 
J=1 /c=l i=l 

i 

771 n m 

A (/ + A, ^ k + tj)u>j(j) k j — C[ 

!c=l j - 1 

i = 1, . . . , m, l 


0 , 

1, • • • ,m 
—u 0 (z(l + t { )), 
l,...,n 


(4.13) 
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and 


?n 

Y2<t>ijL0j\z'(l + tj)\ = 0 , l = 

j = l 

m 

EMAWl = o 5 

3 = 1 


where the intention is that <j>kj ~ + and ~ c^. In this formulation there are actually 

ran + ra + n + 1 equations for the ran + m + n unknowns k = 0, . . - , n, j = 1 , . . . , m and 
c*, k = However, a careful analysis shows that the first 777 equations have a linear 

dependency (the coefficients sum to zero) so that any one of them, e.g., the first, can be 
dropped. This gives a linear system of ran + m + n equations for the unknowns <j>^ 3 and Cfc. 

One need not use the same quadrature rule on the boundary of fi and the cracks. 
Since the solution is smooth on <9fi, we allocate evenly spaced nodes there by t x = i/m, 
i — 0, . . . , ra — 1. The weights are simply u — 1 /ra, corresponding to the trapezoidal rule 
on the closed curve di 7. On each crack <f) will typically have r _ly/2 singularities at the end- 
points and so a quadrature rule is chosen which places more nodes near these singularities 
(but not actually at the endpoints). If the linear crack with endpoint at (a, 6), angle 6 and 
length A is parameterized as (a + t \ cos 0,6-fi iAsin0), 0 < t < 1, then the nodes are chosen 
as (assuming m is even) 


u = /(- 


1/2 


m 


), 


rn 


and t x — 1 - < m -i+ 1 for i = y + 1, . . . , m, where /(x) = 2 q l x q and q is a positive number. 
The parameter q controls the spacing of the nodes with q > 1 causing the nodes to “bunch 
up” near 0 and 1. We have used q — 2.5. The weights are chosen as 


b(ti + ^ 2)1 


t - 1 


U>i = < 


5(<t+i 7 = 2,...,ra- 1 


I 2 ( ^ m — 1 + fm), ^ ^ 


corresponding to a midpoint rule with variably spaced nodes. We have used 60 nodes on 
both d£l and each crack. 

One difficulty with discretizing equations (4.12) is the presence of the logarithmic singu- 
larity in the first kind portion of the equations along the diagonal s = t. We deal with this 
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by using a simple form of product integration based on our quadrature rule. See [2], Section 
3.2, for more details. 


5 The Jacobian and Newton’s Method 


We recall that a central part of our crack reconstruction algorithm is to find a solution of 
the 4 n x 4 n system 


G*(E) = F(£,t/\r*, w^) - F(E*,V> < r fc ,W <rjfc ) 

I 

= F(E,ik,,w„J — f ((/>„„ w,,), i=l,...,n. 


(5.1) 


Here the four components of F(E,^,w) are given by 

r\ f dll)^ 

F(£,V>,^ (,) )= / u(E,0)— — ds i = 1,...,4, 

JdU ov 

where u solves (3.2) and are the functions from (3.9)-(3.1 1). If we use the vector Q £ IR 4n 
to describe the crack configuration E and we parameterize d(Cl \ E) as described in the last 
section, then the discretized version of each of these components becomes 


F{Q, xj>,w) = jr (5.2) 

j = i ov 

where tj is the jth node for Nystom’s method on dtl, t Oj is the corresponding weight. The 
variables Uj are given by Uj = u 0 (tj) + <f>o j, where (f> 0 j form the first part of the solution to 
the system (4.13). 

Our approximate method for the solution of (5.1) is a variation of Newton’s method. For 
that we need an effective method for the calculation of the Jacobian. If q denotes one of the 
components of Q , differentiation of the expression (5.2) yields 


dF_ 

dq 


^(diijdw d ( dw \\ 

E (ity + u % 

J^{d<f>ojdw d ( dw y 




(5.3) 


Here we have assumed that \z'\ is independent of q on dfl and, as mentioned previously, 
we have ignored the functional dependence of the applied current flux rp on E (this in 
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particular gives that ^ = 0). Note that the functions w defined by (3.9)- (3.11) are indeed 
differentiable with respect to the parameters describing the cracks. 

To evaluate from equation (5.3) we therefore need to calculate |^, the derivatives of 
the solution to (4.13) with respect to the parameters which describe the cracks. These can 
be computed with little additional effort. Let the linear system (4.13) be written in the form 

A (Q)MQ) = f(Q) (5.4) 

where A(Q) is the mn + m + n by mn-fm + n matrix appearing in (4.13), f(Q) € JR mn+m+n 
is the right hand side, and <f>(Q) € JR mn+m+ ” is the solution to the system (including the 
constants c). Of course all of these quantities depend on the parameters Q. Differentiation 
of equation (5.4) with respect to any one of the parameters of Q and use of the fact that 
^ = 0 gives 

A ^=f q - d i m=J ^ {Q) ’ (5 - 5) 

so the derivative satisfies a linear equation of exactly the same form as 0 but with a different 
right hand side. Once (5.4) has been solved for <f>(Q) (e.g., by LU decomposition), equation 
(5.5) can be solved for by simply computing the right hand side and reusing the LU 
decomposition. 

Below is a global description of our algorithm for crack reconstruction. We denote by 
E* = {aj}£ =1 the estimated cracks at the ith stage of the algorithm, and by Q* we denote 
the corresponding set of parameters. 

1. Make an initial guess E°, set i = 0. 

2. Select the maximally sensitive two-electrode fluxes 0^ corresponding to the cracks <r£, 
k = 1, . . . , n, in the sense defined in section 3. 

3. Measure (simulate) the boundary voltage data for each flux and compute f(0^,w^), 
k = 1 , . . . , n. 

4. Compute the voltage data for each flux at E = E 2 and compute F(E l , 0^, w^), 
k = 1 , . . . , n. 


21 



5. Compute G^S 1 ) = F(E‘, w <rj: ) — f(VV;p w^), k = 1,. . . ,n. Let G(E) denote the 
vector of all residuals{Gi(E)}’ fc l =1 (a 4n-vector). If G(E’) (= G(Q')) is sufficiently 
small, terminate with the answer E = E*. 

6. Compute the approximate Jacobian J{ = “ DqG(Q)\q = q > v , in the sense described ear- 
lier. 

7. Compute the Newton update SQ ’ by solving J,SQ ‘ = —G(Q’). 

8. Update Q ,+i = Q' + SQ', i = i + 1 , and go to step 2. 

In our implementation of Newton’s method, as in [15], we solve the linear system 

J t SQ' = -G(Q') 

subject to the constraint 

IIW')ll</> 

where A is an appropriately chosen diagonal weighting matrix and p a specified parame- 
ter. This constrained problem is solved in the least squares sense by means of a Lagrange 
multiplier method as outlined in [11]. The constraint on the update markedly improves the 
behavior of Newton’s method far from the solution. 

It should be mentioned that we also constrain the cracks to lie inside the domain ft; 
if the algorithm attempts to move a crack outside ft we perform a simple reflection of the 
endpoint(s) that lie outside to get get the updated crack fully inside ft. The algorithm very 
rarely attempts to move a crack outside the domain (unless the data based on which we 
perform the reconstruction is generated from a collection of cracks, at least which of one lies 
very near the boundary). Our implementation will also allow the cracks to intersect each 
other; the integral equation formulation continues to make sense in this case, although in 
practice one must be careful when dealing with the logarithmic singularities which arise as 
a result of the intersection. We deal with this by again using a form of product integration. 
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6 Computational Experiments 

We have performed a significant amount of computational experimentation with the algo- 
rithm we have developed. In this section we briefly describe five such experiments which 
we find to be representative. The experiments are all performed on simulated data. It is 
our goal to eventually try our algorithm on actual experimental data. A project to build 
an experimental setup and a data gathering device is currently in progress [10]. To go part 
of the way towards reconstruction based on real data, the last of the experiments described 
here pertains to simulated data with a built-in level of noise. 

The domain on which the reconstruction is performed is in all cases the unit disk. The 
graphics by means of which we illustrate the progress of the algorithm is the same for all 
experiments. Each step in the iteration is illustrated by a picture containing two copies of 
the unit disk. The disk on the left depicts the previously estimated locations of the cracks 
(a set of line segments) as well as the boundary locations of the optimally sensitive electrode 
locations. The electrode locations are marked with small circles; the circle corresponding to 
P 0 has been darkened. It is voltage data corresponding to the currents generated by these 
electrode locations that are used for the iterative update. The updated estimates of the 
crack locations are shown as solid line segments in the disk on the right. The true cracks 
which the algorithm seeks to reconstruct, ?.e., the cracks from which the simulated boundary 
voltage data is generated, are shown as dashed line segments (or dashed circular arcs) in the 
disk on the right. 

In our first experiment the simulated data comes from three cracks, two of which are 
very near the boundary. They are each 0.05 units away from the boundary. The cracks have 
lengths 0.25, 0.3 and 0.35. We start the reconstruction procedure by attempting to fit the 
simulated data with data generated from a single crack. As the initial guess we select a crack 
joining the points (-0.2, 0.2) to (-0.2, 0.6). Figure la shows the first step of the iteration. 
After 31 steps the algorithm has found a root for the system of four functionals and reduced 
the residua] error ( | (S Y ( Q ) I ) to 8.94 x 10~ 13 . Figure lb shows step 31 of the iteration. 
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Figure lg: Iteration 19: four cracks fitted to three crack boundary data. 


At this point it is not possible, based on four functionals, to determine whether the 
simulated boundary voltage data comes from just a single crack or several more. The only 
way to determine this is to try to fit more functionals. We take the crack shown in the right 
disk in Figure lb and divide it into two cracks by cutting out a piece 1/10 of the length at 
the center. The two resulting cracks are now used as the initial guess for our algorithm based 
on 8 functionals. The first step of the two-crack iteration is show in Figure lc. After 18 
steps the “two-crack iteration” has found a root to the 8- variable system and has reduced the 
residual error to 6.9 x 10 -16 . The final step is shown in Figure Id. Finally we take the largest 
of the two cracks from the right disk in Figure Id, divide it into two pieces (by cutting out 
the middle 1/10), and give the resulting three cracks as initial guess to our algorithm based 
on 12 functionals. This “three crack iteration” locates the root after 24 steps, reducing the 
residual error to 2.0 x 10 _H . Steps 10 and 24 of this iteration are shown in Figure le and 
Figure If, respectively. If at this point we take and divide one of the three cracks again and 
give the resulting four cracks as an initial guess to our algorithm based on 16 functionals, 
then one of two things is likely to happen: 1) these new four crack will remain essentially 
where the three cracks already are (and the residual will be very small) 2) the algorithm will 
shrink one of the cracks to zero length and the three remaining will stay as before. In the 
latter case the residual will not become quite as small since we impose a lower limit (of 0.01) 
on the length of the cracks. In both cases the behavior clearly indicates that three cracks 
are the right number to fit the data (also for 16 functionals). Figure lg shows step 19 in a 
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“four crack iteration” where the largest crack in Figure If has been divided into two pieces. 
As is evident the latter of the two possibilities from before emerges (the residual after 19 
steps is 5.2 x 10“ 4 ). 

The strategy for gradually increasing the number of cracks as outlined above has emerged 
after a significant amount of numerical experimentation with many sets of simulated data. 
The alternate strategy: to initially guess a sufficient number of cracks and then let the itera- 
tions proceed to convergence has generally shown itself to lead to much slower convergence. 
Our method based on the use of only a few functionals and the adaptive movement of elec- 
trodes has in many of our experiments proven itself to be superior to a least squares fitting 
algorithm (using the entire set of boundary voltage data, but a fixed set of electrodes). Our 
experience with this least squares approach has been that, except for the one crack case, 
it requires a very accurate initial guess in order to converge to the correct solution. For 
multiple cracks the least squares functional appears to contain many local minima. 

This is not to say that our approach may not occasionally be somewhat slower than indi- 
cated by the previous example. We illustrate this with a reconstruction based on simulated 
data from four cracks. Figure 2a shows the first step using a single crack (four function- 
als). After 17 steps our algorithm finds a root with the residual reduced to 1.72 x 10~ 14 . 
However, as seen in Figure 2b (which shows the final iteration) the single crack that is con- 
sistent with the four functionals does not lie near any of the four cracks that were used 
to generate the data. We now divide this single crack into two by cutting out the middle 
1/10. Using the resulting two cracks as initial guess our algorithm based on 8 functionals 
now takes a considerable number of steps before the residual is even reasonably small (and 
the cracks are of reasonable size). Figures 2c and 2d show iterations 99 and 198, respec- 
tively. We take the two cracks after 198 iterations and divide each of them into two using 
the same method as before. The resulting four cracks are provided as initial guess for our 
algorithm based on 16 functionals. Iterations 50 and 115 of this process are shown in Fig- 
ures 2e and 2f, respectively. After 127 iterations a root is found and the residual has been 
reduced to 1.65 x 10 -14 . Even though the algorithm is extremely slow it does ultimately 
converge. It would require an extremely accurate initial guess to get the least squares algo- 
rithm we described before to converge. For comparison figure 2g shows the eleventh and final 
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Figure 2g: Least squares, iteration 11: four cracks fitted to four crack data. 

iterate of a least squares approach. This computation was done using the full boundary 
data and a Levenberg-Marquardt algorithm as outline in [11]. The initial guess used (con- 
sisting of four cracks) is the same as that used for figures 2e and 2f. The least squares 
approach quickly locates a local minimum and terminates. The four cracks that are used 
appear to merge into two. 

Frequently in practice there are either a fairly limited number of well-separated cracks, 
or many cracks clustered in certain locations. Figures 3a and 3b show the first iteration 
and the fifth iteration using our algorithm with one crack (and four functionals) to fit simu- 
lated data coming from 10 cracks located in two clusters. A root is found at the fifth iteration 



Figure 3a: Iteration 1: one crack fitted to ten crack boundary data. 
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guess to the “four crack version” of our algorithm. Figure 4d shows iteration 50 using four 
cracks. The algorithm has not yet located a root. The two longest cracks, which already 
emerge after 13 iterations, provide a reasonable approximation to the curved crack. One of 
the two shorter cracks stays close to the circular arc (and the two large cracks) the other 
one seeks to exit the domain; the program apparently cannot adjust them to find a root, but 
there is a lower bound on their length, so they cannot entirely disappear. 

In the final example we have taken data generated by a two cracks and added 10% 
noise. The noise added is independent and gaussian with a zero mean and standard de- 
viation equal to 1/10 the mean square value of u — Uq on <9f), where u is the potential 




Figure 5a: Iteration 1: two crack fitted to two crack data, 10 percent noise. 



Figure 5b: Iteration 24: two crack fitted to two crack data, 10 percent noise. 
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and Uq is the harmonic function with the same flux as u. In figure 5a we show the first 
iteration. Figure 5b shows the 24th iteration, at which point the algorithm terminates hav- 
ing found a root. Figure 5c shows the result of taking the final crack locations in figure 
5b and using them as an initial guess to a least-squares minimization routine. This rou- 
tine uses a fixed set of electrodes, those from the last iteration of our algorithm. These 
are presumably close to the most sensitive electrode locations for the given cracks. The 
Levenberg-Marquardt routine reduces the total least-squares residual from 0.085 to 0.032 in 
13 iterations and improves the estimate of the crack locations. 



Figure 5c: Least squares, iteration 13: two crack fitted to two crack data, 10 
percent noise. 

7 Summary 

In this paper we have developed a very efficient algorithm for the reconstruction of a collec- 
tion of cracks based on electrostatic boundary measurements. We use a “dual” variational 
formulation for the forward electrostatic problem, and solve this numerically by means of a 
Nystrom’s approximation of the corresponding boundary integral equations. Our reconstruc- 
tion is based on adaptively changing the current patterns, so as to maximize the sensitivity 
of the measured voltage differences. Our reconstruction is based on a limited set of aver- 
ages of the boundary voltage measurements as opposed to the entire set of measurements; 
this should lead to greater efficiency and less rigidity. The algorithm is currently entirely 
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two-dimensional, but it should be very interesting to extend it to three dimensions. At this 
point we have only investigated the behavior of the algorithm when used on “synthetic” 
data (including data with noise). It should be very interesting to apply the algorithm to 
data coming from “real” experiments. Frequently cracks appear as clusters of many small 
(microscopic) cracks; our algorithm, when applied to data generated by clusters of small 
cracks, often very successfully locates a set consisting of a few, well-separated cracks. In this 
context it should be extremely interesting to analyze in what sense this reflects the behavior 
of the forward problem. To be more specific: it should be interesting to study in what sense 
a cluster of small (microscopic) cracks in an appropriate limit approaches a single (lumped) 
macroscopic crack. 
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